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Abstract 

We consider a model where a population of diffusively coupled limit-cycle oscillators, 
described by the complex Ginzburg-Landau equation, interacts nonlocally via an 
inertial field. For sufficiently high intensity of nonlocal inertial coupling, the system 
exhibits birhythmicity with two oscillation modes at largely different frequencies. 
Stability of uniform oscillations in the birhythmic region is analyzed by means of 
the phase dynamics approximation. Numerical simulations show that, depending on 
its parameters, the system has irregular intermittent regimes with local bursts of 
synchronization or desynchronization. 



1 Introduction 

Beginning with the pioneering contributions by Kuramoto [1] and Winfree 
[2] , studies of synchronization and spatiotemporal chaos (turbulence) in pop- 
ulations of coupled active oscillators have developed into a classical field of 
research. In chemical systems, each oscillator represents a reaction element 
and coupling between such elements is usually due to diffusion of reactants 
in the system. This coupling extends only within a short diffusion length and 
is therefore local. The complex Ginzburg-Landau equation is the canonical 
model for oscillatory systems with local coupling near a supercritical Hopf 
bifurcation. In surface chemical reactions [3], reaction elements are however 
additionally interacting through the gas phase where instantaneous complete 
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mixing occurs. As a result, global coupling between chemical oscillators arises 
[4]. Moreover, global delayed feedbacks through the gas phase could be artifi- 
cially introduced to control turbulence in surface reactions [17,18,19,26,27,28]. 
A special class of systems are arrays of active oscillators that do not directly 
interact one with another, but are all coupled to a single diffusing field (in 
biochemistry, an example of such a system would be provided by arrays of al- 
losteric enzymes interacting through a chemical messenger [5]). Adiabatically 
eliminating this field, models with nonlocal coupling between the oscillators, 
described by a finite-range integral, were derived [6,9]. Investigations of non- 
locally coupled oscillator arrays have revealed that such form of coupling does 
not always synchronize neighbouring oscillators and discontinuous distribu- 
tions with scaleless fractal structure can develop [8,9]. Spatiotemporal chaos 
can develop in such systems even in the Benjamin-Feir stable regime [12] and 
spiral waves with phase-randomized cores can exist here [15,14]. 

Furthermore, arrays with both local and nonlocal coupling between oscillators 
are possible. This situation is characteristic, for example, for surface chemical 
reactions. In such systems, diffusion provides local coupling between neigh- 
bouring surface oscillators, whereas much faster heat conduction is responsible 
for nonlocal coupling between them. 

Recently, Tanaka and Kuramoto [16] have shown how, in the vicinity of a 
supercritical Hopf bifurcation, the description of arrays of nonlocally coupled 
oscillators can be reduced to the complex Ginzburg-Landau equation with 
nonlocal coupling. Because of the critical slowing down near the bifurcation 
point, the coupling is effectively instantaneous in the reduced description. 

In realistic systems, which are not too close to the Hopf bifurcation, one can 
however also encounter the opposite situation, with a very slow inertial field 
giving rise to nonlocal coupling. The aim of our study is to analyze what 
new effects, primarily due to slow nonlocal coupling, are possible in this class 
of systems. For our investigations, we have chosen an abstract model of the 
complex Ginzburg-Landau equation interacting with an additional slow field, 
so that the coupling is nonlocal both with respect to space and time. We 
expect that the behaviour found in this general model would be typical for a 
broad class of realistic systems. 

The principal effect of coupling inertiality is that birhythmicity, leading to 
chaotic intermittency, can develop in such systems. In contrast to the birhyth- 
micity in reaction-diffusion systems near a pitchfork-Hopf bifurcation (see 
[20,21]), the two oscillatory states are characterized here by very different 
frequencies and oscillation amplitudes. For the rapid oscillatory state, the ad- 
ditional field responsible for the inertial long-range coupling is almost absent 
and only diffusive local coupling between the oscillators is important. In slow 
oscillations, the elements are however entrained by the long-range inertial field. 
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Depending on the parameters, two new kinds of intermittency are found here. 
When rapid oscillations are modulationally unstable and give rise to turbu- 
lence, spatial regions occupied by slow entrained oscillations spontaneously 
develop and die out in the medium. They can be interpreted as bursts of syn- 
chronization on a turbulent background. On the other hand, relatively small, 
irregularly evolving islands filled with rapid turbulence can persist on the back- 
ground of slow almost uniform oscillations. They can be therefore described 
as bursts of desynchronization on the background of regular slow oscillations. 
Our analysis of such phenomena is based on the phase dynamics approxima- 
tion, derived for the birhythmic system. It is complemented by ID and 2D 
numerical simulations. 



2 The Model 



The investigated model describes a system of diffusively coupled active oscil- 
lators coupled to an additional inertial diffusive field. Introducing the local 
complex oscillation amplitude rj(x,t) and denoting as z{x,t) the additional 
complex- valued diffusive field, we have therefore a system of two equations 

77 = (1 + iuj)r] - (1 + ia)\r]\ 2 r] + (1 + i/3)V 2 ry + K(z - rj) (la) 
rz = rj-z + l 2 V 2 z (lb) 

The additional last term in the complex Ginzburg-Landau equation (la) 
takes into account coupling of the oscillatory subsystem with the field z whose 
evolution obeys equation (lb); K is the respective coupling constant. The 
equations are brought into a dimensionless form by choosing the characteristic 
diffusion length in the oscillatory subsystem as the length unit and taking the 
characteristic relaxation time scale of the oscillators as the time unit. The 
parameters r and I determine characteristic time and length scales of the 
additional field z. We assume that this field is inertial (r 3> 1) and slowly 
varying in space (/ 3> 1). 

The linear equation (lb) can easily be solved as 

/oo rt 
/ G(x -x',t- t f )r](x\ t')dx'dt'. (2) 
-00 Jo 

where the kernel is given by 

G(M)= 5^ exp ("is"9- < 3 > 

Substituting this into equation (la), we obtain an equivalent integro-differential 
form of the considered model, 
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?7 = (l + iu;)?7- (1 + ia)|?7| 2 ?7 + (1 + i/9)V 2 ?7 

/oo rt 
/ G{x -x',t- t f ) [r)(x\ t f ) - 7](x, t)] dx f dt f (4) 
-oo JO 

We see that, besides of diffusion, the model also includes an additional cou- 
pling, nonlocal both with respect to space and time. 

Note that very close to a supercritical Hopf bifurcation all processes become 
fast as compared to the relaxation time scale of individual oscillators, because 
this time is inversely proportional to the distance from the bifurcation point 
(critical slowing down). However, it is known that the complex Ginzburg- 
Landau equation yields qualitatively correct description even relatively far 
from the bifurcation. Equations (1) and (4) can be viewed as providing a 
simple model of an oscillating system coupled to an inertial diffusive field. 



3 Birhythmicity 



The system described by equations (1) is birhythmic, i.e. it can have two 
different oscillatory uniform states. Assuming that rj(t) = pe~ llt and z(t) = 
re~ llt and substituting this into equations (la) and (lb), we obtain a cubic 
equation for the oscillation frequency 7, 

r 2 7 3 + t 2 (uo - a + aK)rf + (1 + tK)j + uu-a = 0. (5) 



When the frequency 7 is known, the oscillation amplitude p of the field 77 is 
given by 

1 + r 2 7 2 

The respective oscillation amplitude of the field z is 



1 - TT-h- <6) 



Depending on its parameters, equation (5) can have either one or three real 
roots. The three roots correspond to three possible modes of uniform oscilla- 
tions with frequencies 71,2,3, such that 71 < 72 < 73. It can be checked that 
oscillations with the middle frequency 72 are always unstable. In contrast to 
this, uniform oscillations with frequencies 71 and 73 are possible (but may 
still be unstable with respect to nonuniform perturbations, see the discussion 
below) . 

Figure 1(a) shows the region in the parameter planes (a — K) and (r, K) 
where birhythmicity exists. The birhythmicity is possible only for sufficiently 
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Figure 1. Birhythimicity regions (gray) for the model (1) in the parameter planes 
(a — uj,K) and (t,K). The fixed parameters are r — 10 in the first diagram and 
uo — 2, a — 2.5 in the second diagram. The boundaries of the displayed regions do 
not depend on the parameter (3 of the model. 

strong coupling K. Moreover, it develops only if the characteristic time r of 
the field z is sufficiently large (Fig. 1(b)). 

Numerical solutions of equation (5) are displayed in Fig. 2. If the param- 
eter a is kept constant and the coupling strength K is varied (Fig. 2(a)), 
birhythmicity is found inside the interval 0.4 < K < 1. When K > 1, 
the system has a stable stationary state coexisting with oscillations. The de- 
pendence of the effective oscillation frequency 7 on the difference a — u for 
K = 0.5 is shown in Fig. 2(b). Two stable limit cycles coexist within the in- 
terval -0.1 < a- uj < 0.8. 

The difference between the two limit cycles becomes clear if we compare the 
amplitudes of oscillations of the fields rj and z which are presented below in 
the same figures. As we have already observed, at K = only one oscillating 
mode is present, as we are in this limit reduced to the standard CGLE. Thus, 
we can expect that the branch starting at K = (which is the upper one in 
the bifurcation diagrams of Fig. 2) would have the strongest similarities to 
the CGLE, even though for finite K it would be affected by the coupling to 
the second field. This branch is characterized by a large amplitude of the field 
r] ) which never gets significantly smaller than 1. The corresponding frequency 
is quite small, starting from 7 = 0.5 at K = and slowly decreasing as K 
increases. Below, we refer to this attractor as corresponding to the slow limit 
cycle. 

In contrast to this, oscillations in the lower branch have a smaller amplitude 
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Figure 2. Bifurcation diagrams. Frequencies 7 of uniform oscillations, together with 
the respective amplitudes p = \rj\ and r — \z\ as functions of K, a — u and r for 
slow (bold line), rapid (dashed) and absolutely unstable (dotted) uniform oscillation 
modes. The parameters are (a) u = 2, a = 2.5, r = 10 (b) u = 2, K = 0.5, r 10, 
(c) cj = 2, a = 2.5, if = 0.5. 

p, which decreases for higher coupling strengths K and eventually vanishes at 
K — 1. Note that 7 is negative for this mode and therefore such oscillations 
are counter-rotating with respect to the upper branch. The frequency varies 
significantly for this mode. Inside the birhythmicity interval, the frequency 
of such rapid limit cycle is always much larger (in its magnitude) than the 
frequency of the slow mode. A characteristic feature of the rapid limit cycle 
is that the amplitude r of the field z is very small here. This is because the 
oscillations of the field r\ are so rapid here that the inertial field z cannot 
respond to them. Hence, the coupling term Kz in the system (1) is almost 
negligible for the rapid limit cycle so that the nonlocal coupling is ineffective 
and oscillators interact only via the local diffusional coupling in this mode. 



4 The Phase Dynamics Approximation 

Above, uniform oscillations in the model have been discussed. Now, we want to 
consider the properties of almost uniform oscillations, characterized by small 
phase gradients. Because the system is invariant with respect to the uniform 
shift of all phases, the evolution of phase distributions with small phase gra- 
dients should be slow, in contrast to fast local relaxation of the amplitude 
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perturbations. Therefore, a reduced description, known as the phase dynam- 
ics approximation, can be constructed in such cases [1]. When birhythmicity 
is present, the coefficients of the phase dynamics equation are different for the 
two oscillation modes. 



The derivation of the phase dynamics equation for this system is lengthy and 
has been performed by using a computer program for algebraic calculations 
(Mathematica) . Below in this section we show the scheme of derivation and 
discuss the results. The complete derivation is presented in Appendix A. 



We define the phases (cp and p) and amplitudes (p and r) of complex fields 
r] and z as 77 = pe 1 ^, z = re lip . It is convenient to use instead of (p and p the 
variables ip = (p — p and 6 = (p + p. Note that the phase sum 6 is the slow 
variable in this system. The phases of the two fields are rigidly correlated in 
the uniform state and their difference ip undergoes rapid relaxation, similar to 
oscillation amplitudes p and z. 



In terms of the new variables, the model (1) takes the form 



p = p — p 3 + Kr cos ip — Kp + p 

P T I 2 

r = - cos ip — r VB 2 
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Suppose now that the birhythmic system is close to one of its two oscillatory 
states, characterized by some amplitudes p and r , and the relative phase 
^o, whereas the phase sum 6 remains arbitrary. If the phase sum is slowly 
varying in space, the variables p, r and ip would only slightly deviate from their 
equilibrium values. Therefore, we can introduce small perturbations of such 
variables, p = p + 5p, r = r + 5r, ip = ip + Sip, and linearize the equations 
for with respect to such perturbations. The linearized evolution equation for 
the phase 9 is 
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(9) 



The linearized equations for the variables 5p, 5r, Sip can be solved in the adia- 
batic approximation, because these fast variables are adjusting to slow varia- 
tion of 6. Substituting the result into (9), a closed evolution equation for the 
phase variable 6 is derived. This equation has the form: 



The analytical expressions for the coefficients Co, C\ are given in Appendix 
A. Note that inside the birhythmic region there are two uniform oscillation 
modes with different p , ^o, and ip , and thus with different values of these 
coefficients. 

The most important coefficient is C 2 because its sign controls stability of uni- 
form oscillations with respect to phase modulation (the Benjamin-Feir insta- 
bility). If C2 > 0, uniform oscillations are stable, if C2 < they are unstable. 
We have computed this coefficient as functions of several model parameters, 
by using the obtained analytical expressions. The computed dependences are 
shown in Fig. 3. Two different branches correspond to the two different limit 
cycles. The upper curves in Fig. 3(a,b,d) are for the slow mode, whereas the 
lower curves are for the rapid oscillations. In Fig. 3(c), the upper curves cor- 
respond to the slow mode for larger values of the parameter /. 

Examining these plots, several observations can be made: Slow mode oscilla- 
tions are always stable inside the birhythmic region, even when the system 
without coupling is Benjamin-Feir unstable (i.e., 1 + a(3 < and therefore 
C2 < at K = 0). In contrast to this, rapid oscillations are always unstable 
when the condition 1 + a (3 < is satisfied. Moreover, they are unstable near 
the boundary of the birhythmic region, at which they first appear, even when 
this condition is violated. Increasing the diffusion length / of the field z favors 
instability of rapid uniform oscillations. 



e = c + c f 1 (ve) 2 + c 2 v 2 e 



(10) 
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Figure 3. The coefficient C2 of the phase dynamics equations as functions of K, a, / 
and r for several different values of the parameter (3. The upper branches in (a,b,d) 
and for the larger values of / in (c) correspond to slow oscillations, the lower branches 
are for the rapid oscillations. The oscillations are unstable if C2 < 0. The parameters 
are (a) a = 2.5, I = 10, r = 10, (b) K = 0.5,/ = 10, r = 10, (c) a = 2.5, K = 0.5, 
t 10, (d) a 2.5, K 0.5, I 10; for all curves u = 2. 

5 Numerical Simulations 

To investigate nonlinear spatiotemporal dynamics in the considered model, 
simulations of equations (la) and (lb) have been performed. For numerical 
integration of these equations, the fourth-order Runge-Kutta algorithm has 
been used. The mesh size for space discretization and the time step have been 
chosen to optimize the computational time for each parameter choice. Since 
the diffusion length and the diffusion constant of the oscillatory field have 
both been chosen to be equal to unity, Ax varied between 0.3 and 0.5, while 
At could vary between (Ax) 2 /5 and (Ax) 2 /2. Both one- and two-dimensional 
systems were investigated. No-flux boundary conditions were employed. Initial 
conditions varied depending on a particular simulation. 

To display the results of one-dimensional simulations, space-time diagrams 
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(a) Re(7 ? ) (b) \ V \ 

Figure 4. Spatiotemporal diagrams showing front propagation in the birhythmic 
system with uj = 2, a = 2.5, (3 = 1, K = 0.5, / 10, and r 10. The total system 
length is L = 120 and the total displayed time interval is T = 300. Bright regions 
correspond to small values of displayed variables. Time runs from left to right along 
the horizontal axis, the spatial coordinate is varied in vertical direction. 

have been constructed. In such diagrams, time always runs along the horizon- 
tal axis and the vertical axis is corresponding to the spatial coordinate. In 
two dimensions, snapshots at some time moments and space-time diagrams 
showing evolution of the pattern along some linear cross-section are presented. 
The patterns are shown by using a gray scale where the white color encodes 
the smallest and the black color the largest values of the displayed variable. 

First, the behavior of fronts separating spatial regions with two different os- 
cillatory states in the birhythmic system has been studied (Fig. 4). The front 
travels towards the slow oscillating region, in such a way that the system ends 
up with uniform oscillations of the rapid mode. Its propagation is character- 
ized by periodic appearance of amplitude defects which occur when the phase 
difference between the two oscillating regions equals 2tt (see Fig 4(b)). The 
front moves at a constant velocity which depends on the parameters of the 
system. Fig. 5(a) shows the dependence of the front velocity on the parameter 
a for several values of the coupling strength K. In Fig. 5(b), the dependence 
of the front propagation velocity on the frequency difference A7 = 73 — 71 
of two uniform oscillatory modes is presented. The front always moves faster 
when this difference is increased. 

As follows from the stability analysis of uniform oscillatory states (see, e.g., 
Fig. 3(b)), rapid oscillations should be unstable with respect to phase modu- 
lation near the birhythmicity boundary and complex spatiotemporal regimes 
can be expected there. We have performed a series of numerical simulations, 
exploring pattern formation in the system within the intervals 2.55 < a < 2.8 
and —1 < (3 < 2.5, while keeping fixed the other parameters (u = 2,K = 
0.5, 1 = 10, and r = 10). Each simulation started from independently chosen, 
completely random initial conditions. 

Figure 6 is a schematic summary of these numerical investigations. The solid 
line is the stability boundary of rapid oscillations (C2 = 0); rapid oscillations 
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Figure 5. Front velocity as functions of the parameter a and difference A7 of oscil- 
lation frequencies in the two modes in the birhythmic regime. The fixed parameters 
are u = 2, (3 = 1, 1 = 10, and r — 10. 




Figure 6. Schematic phase diagram. Numerical simulations of the one-dimensional 
model were performed at the values of a and (3 indicated by symbols in this dia- 
gram. Depending on the observed properties of patterns, the diagram is divided into 
regions 1 to 6. The circles indicate the values of these parameters used to produce 
a typical patterns for the corresponding region, displayed in Fig. 7. The solid line 
shows the stability boundary of rapid uniform oscillations, given by the condition 
C2 — 0. Other parameters are u — 2, K — 0.5, / — 10, and r — 10. 

are unstable below this boundary. Slow oscillations are always stable in the 
considered region. Simulations have been performed for a set of parameter 
values indicated by symbols in this diagram. 

Comparing properties of observed patterns, we can qualitatively divide them 
into several groups occupying different regions in the parameter plane. Each 
group is marked by its own symbol and the respective regions in Fig. 6 are 
numbered. Examples of typical spatiotemporal patterns, observed in one- 
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dimensional simulations within each region, are displayed in Fig. 7. These 
simulations were all performed for a system of the total length L = 200 and 
for a total time varying between 500 to 3000 time units. To clearly present the 
observed patterns, only parts of their entire evolution are however displayed. 

In region 1, the slow limit cycle is Benjamin- Feir stable, while the rapid limit 
cycle is unstable. The system is found here in the regime of fully developed 
amplitude turbulence characterized by creation of multiple amplitude defects. 
Increasing parameter /?, we enter a region 2 closer to the stability boundary 
where turbulence becomes intermittent. We observed the emergence of larger 
groups of synchronized oscillators which become able to reach the rapid limit 
cycle and perform harmonical oscillations for several periods. Inside the small 
parameter region 3 near the stability boundary, the frequency and the ampli- 
tude of oscillations correspond almost everywhere to those of the rapid limit 
cycle. However, the system does not undergo complete synchronization. Lines 
of amplitude defects travelling through the system act as wave sources here. 

In regions 5 and 7, both lying above the stability boundary of rapid oscil- 
lations, uniform oscillations are observed starting from random initial condi- 
tions. The oscillations are rapid inside region 5 and slow inside region 7. These 
two domains are separated by regions 4 and 6 in the parameter plane, where 
competition between two oscillation mode takes place leading to complex spa- 
tiotemporal regimes. 

The patterns inside region 4 are characterized by a background of rapid chaotic 
oscillations, with small amplitude and numerous amplitude defects. On this 
highly desynchronized background, bursts of synchronization emerge. Such 
bursts consist of large groups of elements which suddenly start to oscillate 
together with a large amplitude and a small frequency corresponding to the 
stable slow limit cycle. However, these groups do not keep synchronized over 
a long time: after less than one oscillation period the turbulence overwhelms 
again. As already noticed in our discussion of birhythmicity, the amplitude of 
the coupling field z is very small in the rapid limit cycle, while in the slow 
limit cycle it gets a larger value comparable to the amplitude of rj. Analyzing 
patterns in region 4, it can be seen that during the synchronization bursts, the 
amplitude \z\ is indeed much larger than in the turbulent background. Note 
that such patterns with synchronization bursts have also been observed for 
some parameter values below the stability boundary of rapid uniform oscilla- 
tions. 

Region 6 lies near region 7, where slow uniform oscillations are winning the 
competition. Here, the patterns can be described as exhibiting bursts of desyn- 
chronization on the background of slow uniform oscillations. Inside such bursts, 
the coupling field z is strongly decreased in amplitude and only short-range 
diffusive coupling between oscillators is effective. 
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Figure 7. Spatiotemporal diagrams displaying evolution of Re(?7), p — \rj\ and r = \z\ 
in typical patterns observed in the regions 1-6 for the one-dimensional system. The 
respective parameter values are given in Fig. 6. The displayed space and time in- 
tervals are (1) L = 100, T = 166, (2) L = 100, T = 166, (3) L = 100, T = 250, (4) 
L = 200, T = 500, (5) L = 100, T = 250, (6) L = 200, T = 500, and (7) L = 200, 
T = 250. The contrast level is adjusted individually in each plot to ensure best 
visualization of the pattern. 
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Figure 8. Selected snapshots of the phase portraits for the pattern 6 in Fig. 7. The 
circles show two different limits cycles in the birhythmic system. 

To further illustrate this regime, phase portraits of the patterns in region 6 
were constructed at different time moments. Fig. 8 shows several subsequent 
snapshots from the recorded video. Each point in a phase portrait represents 
a state of one of the oscillators. The spatial information (i.e., spatial distances 
between the displayed oscillators) is lost in this representation, but the motions 
performed by the oscillator population are seen more clearly. The two circles 
indicate the two coexisting limit cycles. The cycle with a smaller oscillation 
amplitude is rapid, whereas the cycle with a larger amplitude is slow. 

We can see from these snapshots that the system is usually divided into two 
groups of oscillators, occupying the two coexisting attractors (limit cycles). 
Some oscillators are found in the vicinity of the origin rj = 0. They correspond 
to amplitude defects generated at the border between spatial regions with 
different oscillation frequencies (bright dots in the spatiotemporal plot of \rj\ 
for pattern 6 in Fig. 7). Such defects are not present at all times, and the 
respective points are not, for example, seen in the phase portraits at t — 16.7 
and t = 281.7. The oscillators sitting on the rapid limit cycle with the smaller 
amplitude belong to the desynchronization bursts. 

For selected parameter values, two-dimensional simulations of the system 
have additionally been performed. Figure 9 displays the behavior of a two- 
dimensional system at the same parameters as for the pattern 3 in Fig. 7. The 
space-time diagrams in the right panels show the pattern development along a 
horizontal cross section. Starting from random initial conditions, a population 
of rotating spiral waves develops. After a transient, a stable configuration of 
rotating spirals is reached in two dimensions. 



14 




Figure 9. Multiple spiral waves in a two-dimensional system of size 120 x 120; 
the same parameter values as for the pattern 3 in Fig. 7. The right panels are 
space-time plots showing evolution of the pattern along one horizontal cross section; 
the displayed time interval is T — 480 

Figure 10 shows two-dimensional patterns corresponding to synchronization 
bursts (pattern 4 in Fig. 7). Inside such a burst, occupying for example the left 
central region in the left panels, the coupling field z is increased in magnitude. 
The space-time diagrams (right panels) reveal that the synchronized regions 
with slow oscillations have only relatively short lifetimes and are replaced 
by irregular rapid oscillations. However, they appear again and again in the 
course of time. 

A behavior, which can be described as bursts of desynchronization, was found 
in our two-dimensional simulations even outside of the birhythmicity region, 
where only slow uniform oscillations are possible (Fig. 11 and 12). We started 
here with the initial condition, which is commonly used to generate rotating 
spirals. A rotating spiral was indeed first formed. However, some instability 
has then developed beginning from the its central region, where the oscillation 
amplitude z was decreased (see Fig. 11). The development of this instability 
has resulted (Fig. 12) in complete destruction of the spiral and the appearance 
of relatively small chaotic domains on the background of almost uniform slow 
oscillations. 

Typical snapshots of spatial distributions of the variables Re(r/), and \z\ in 
this pattern are displayed in the left panels in Fig. 12. Inside the domains, the 
coupling field z is reduced in magnitude and these small spatial regions are 
filled with irregular rapid variations of the complex field rj. Thus, they can be 
classified as desynchronization bursts. The space-time diagrams through one 
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Figure 10. Bursts of synchronization in a two-dimensional system of size 120 x 120; 
the same parameter values as for the pattern 4 in Fig. 7. The right panels are 
space-time plots showing evolution of the pattern along one horizontal cross-section; 
T = 480. 




Figure 11. Instability of a rotating spiral. Simulation for a two-dimensional system 
of size 120 x 120 with parameters u = 2, a = 3, (3 = 1, K = 0.4, / 10, and r 10. 
Subsequent snapshots of the field Re(r/) at times T = 2.8, 10.0, 15.2, 37.2, 96.8, 
119.6, 130.0, and 139.2 are presented. 

horizontal cross section of the medium (right panels in Fig. 12) show that such 
domains can spread or shrink, and travel through the medium. The edges of 
these structures are marked by the appearance of amplitude defects where \rj\ 
is close to zero. 

Remarkably, if we keep fixed all other parameter values, but eliminate diffu- 
sional coupling between the oscillators (the Laplacian term in equation (la)), 
evolution starting from the same initial conditions leads to spiral waves with 
phase-randomized cores, similar to those which were previously considered 
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Figure 12. Bursts of desynchronization in a two-dimensional system of size 
120 x 120; continuation of the simulation presented in Fig. 11. The parameters 
are u = 2, a = 3, = 1, K = 0.4, 1 = 10, and r = 10; T = 720 

[15,14]. Thus, the inclusion of diffusive coupling has a strong effect on pat- 
tern formation in the system. Though patterns resembling spiral waves with 
phase-randomized core are indeed initially developing in the diffusively cou- 
pled case, they are unstable and, after a transient, lead to the development of 
intermittent spatiotemporal regimes with desynchronization bursts. 



6 Discussion 

Our study can be viewed as an extension of a series of detailed investiga- 
tions of systems with nonlocal coupling between oscillators performed by Y. 
Kuramoto with his coworkers. We have chosen the model (1), which belongs 
to the previously discussed class, and focused our analysis on the effects of 
coupling inertiality and diffusion in this model. The principal effect of the in- 
ertiality in nonlocal coupling is that the system becomes birhythmic and has 
uniform slow and rapid oscillations as two different coexisting attractors. 

For rapid oscillations, the inertial field z responsible for nonlocal coupling is 
not excited and the system is essentially reduced to an array of only diffusively 
coupled oscillators with amplitudes rj. In contrast to this, the nonlocal coupling 
field z is involved in slow oscillations and nonlocal coupling between oscillators 
is in operation under these conditions. 

Our stability analysis based on the phase dynamics equation has shown that, 
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though slow uniform oscillations are always stable in the considered model, 
rapid uniform oscillations may become modulationally unstable near the birhyth- 
micity boundary. Traveling fronts, separating spatial regions with two different 
oscillation modes, have been seen in our numerical simulations. Note that the 
fronts separating oscillating and non-oscillating regions have been recently 
observed in the vicinity of a subcritical Hopf bifurcation [24] . 

Irregular intermittent spatiotemporal regimes of two different kinds have been 
numerically observed. In the patterns with synchronization bursts, the back- 
ground is occupied by rapid chaotic oscillations. On this background, short- 
living small domains with slow and rather uniform oscillations spontaneously 
develop. Inside such domains, the nonlocal coupling field is much larger in mag- 
nitude than for the rapidly oscillating background. The characteristic size of 
such domains was comparable to the nonlocality radius in the model. This kind 
of intermittent turbulence is to some extent similar to the behavior observed 
for a system near the Andronov homoclinic bifurcation, exhibiting bistability 
between a stationary and an oscillatory state [22,23]. However, the developing 
domains are filled here not by a stationary state, but by slow oscillations. 

In the spatiotemporal patterns with desynchronization bursts, the background 
is occupied by slow oscillations which are almost uniform. On this background, 
a cascade of desynchronization bursts develops. Each burst represents a local- 
ized turbulent spot filled with rapid irregular oscillations. Inside it, the non- 
local coupling field is reduced in magnitude. Such patterns correspond to the 
regimes of intermittent localized turbulence previously seen for the complex 
Ginzburg-Landau equation with global coupling [18], in the experiments with 
global delayed feedback in surface chemical reactions [26] and in the respective 
realistic simulations [27]. In some cases, they develop in the situations when 
rotating spiral waves with phase-randomized cores have been seen [15,14] in 
systems where only nonlocal coupling was present. Birhythmicity as a cause 
of weak turbulence has recently been investigated in a biological model of 
glycolitic oscillations [25]. 

Stimulating discussions with Prof. Y. Kuramoto are gratefully acknowledged. 
We thank M. Stich for the discussion of questions related to birhythmicity. 



Appendix A 



Starting from system (8) we linearize the first three equations around the 
values po, To, of the uniform oscillations. The system that we obtain can be 
written as 
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Sp = a 1 Sp + hSr + aStf) + c?iV6 2 + eiV 2 6 (11a) 

Sr = a 2 Sp + b 2 Sr + c 2 8ip + d 2 V6 2 (lib) 

Sip = a 3 5p + b 3 Sr + c 3 5ip + d 3 V6 2 + e 3 V 2 6 (11c) 

= a 4 Sp + b 4 Sr + c 4 5^p + d 4 V6 2 + e 4 V 2 6 + f 4 (lid) 

where the coefficients are given in the following table 



ai = 1 - 3pl - K 


a 3 = 2«p + 




sin(^o) 


b\ = K cos (ibn ) 


6 3 = sin(^o) 


Po rr^. 




Cl = -Kr sin(^o) 


c 3 = cos(^o) 


— J(IQ- _L_ _£o_ 

Po rr _ 




di = -? 


ds = -l 2r f T 


ei = §Po/5 


^3 = K 1 -?) 


a 2 = I cos(^o) 


a 4 2apo + 


L Pn r ^o 


sin(^o) 




64 = sin(^o) 


_K _po_~ 

Po rr 2 _ 




c 2 = -^sin(^o) 


c 4 = cos(^o) 


_j^m 1 _£q_ 

Po rr _ 




^2 = -J 2 £ 


d 4 = f 


e 2 = 




f 4 = -oj + apl + [-K% + sin(^o) 



Now, since we are in the approximation where p, r, ^ adjust adiabatically to 
0. we can assume Sp = Sr = Sip = 0, so that we get from (lid) 



^ b 3 c 2 d x - b 2 c 3 d 1 - b 3 c x d 2 + hc 3 d 2 + &2Cirf 3 &ic 2 4 ^q2 
-a 3 b 2 c x + a 2 6 3 ci + a 3 b x c 2 - a r b 3 c 2 + a 2 6ic 3 - a x b 2 c 3 

_! & 3 c 2 ei - b 2 c 3 e 1 + 6 2 Cie 3 - b x c 2 e 3 ^ 2q 

-a 3 6 2 ci + a 2 6 3 ci + a 3 &ic 2 - ai& 3 c 2 + a 2 b r c 3 - a x b 2 c 3 

^ _ a 3 c 2 d x - a 2 c 3 d x - a 3 c x d 2 + a\C 3 d 2 + a 2 c Y d 3 - a\C 2 d 3 ^^ 2 

a 3 6 2 ci - a 2 b 3 ci - a 3 b r c 2 + a r b 3 c 2 + a 2 b\C 3 - a r b 2 c 3 
_^ a 3 c 2 ei - a 2 c 3 ei + a 2 c x e 3 - a x c 2 e 3 ^ 2q 

a 3 6 2 ci - a 2 b 3 ci - a 3 b x c 2 + a r b 3 c 2 + a 2 &ic 3 - ai& 2 c 3 
— a 3 Wi - « 2 6 3 (ii - a 3 b r d 2 + ai& 3 d 2 + a 2 b r d 3 - a r b 2 d 3 ^^ 2 
-a 3 b 2 ci + a 2 6 3 ci + a 3 &ic 2 - a x b 3 c 2 - a 2 b x c 3 + ai& 2 c 3 

^ a 3 & 2 ei - a 2 6 3 ei + a 2 b x e 3 - a x b 2 e 3 ^ 2q 

-a 3 b 2 ci + a 2 b 3 ci + a 3 b r c 2 - a r b 3 c 2 - a 2 b r c 3 + a x b 2 c 3 

These expressions can be put into the equation for to get 
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e = c + c 1 (ve) 2 + c 2 v 2 e 

where 



(15) 



CW4 (16) 
Ci = rf 4 + [c 4 (a 3 6 2 rfi - ^263^1 - ^361^2 + CLib 3 d 2 + a 2 b x d 3 - aib 2 d 3 ) 
-b 4 (a 3 c 2 di - a 2 c 3 di - a 3 c x d 2 + a x c 3 d 2 + a 2 Cid 3 - aic 2 d 3 ) 
+a 4 (b 3 c 2 d 1 - b 2 c 3 d x - b 3 c ± d 2 + b x c 3 d 2 + b 2 c ± d 3 - 6iC 2 rf 3 )] / 
{-a 3 b 2 ci + a 2 b 3 ci + a 3 b r c 2 - a x b 3 c 2 - a 2 b x c 3 + a r b 2 c 3 ) (17) 
C 2 = e 4 + [c 4 (a 3 6 2 ei - a 2 b 3 e x + a 2 6ie 3 - a r b 2 e 3 ) 
-6 4 (a 3 c 2 ei - a 2 c 3 ei + a 2 c r e 3 - aic 2 e 3 ) 
+a 4 (6 3 c 2 ei - b 2 c 3 e r + b 2 c r e 3 - bic 2 e 3 )] / 

{-a 3 b 2 ci + a 2 b 3 ci + a 3 b t c 2 - a x b 3 c 2 - a 2 b x c 3 + aib 2 c 3 ) (18) 
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